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This paper addresses the Bayesian calibration of dynamic models with parametric and struc- 
tural uncertainties, in particular where the uncertain parameters are unknown/poorly known 
spatio-temporally varying subsystem models. Independent stationary Gaussian processes with 
uncertain hyper-parameters describe uncertainties of the model structure and parameters 
while Karhunnen-Loeve expansion is adopted to spectrally represent these Gaussian pro- 
cesses. The Karhunnen-Loeve expansion of a prior Gaussian process is projected on a gener- 
alized Polynomial Chaos basis, whereas intrusive Galerkin projection is utilized to calculate 
the associated coefficients of the simulator output. Bayesian inference is used to update the 
prior probability distribution of the generalized Polynomial Chaos basis, which along with 
the chaos expansion coefficients represent the posterior probability distribution. Parameters 
of the posterior distribution arc identified that quantify credibility of the simulator model. 
The proposed method is demonstrated for calibration of a simulator of quasi-one-dimensional 
flow through a divergent nozzle. 
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Figure 1. Conceptual Architecture of Bayesian Calibration for Uncertain Models 



1. Introduction 



With the present ubiquitous use of computer simulators for the scientific inves- 
tigations, uncertainty quantification and calibration of the simulator models is 
identified as an important area of research [H-0]. Significant developments of the 
last decade have established the Bayesian framework as a preferred method for 
uncertainty quantification and calibration of computer simulators [cj-tl^]. This pa- 
per explores a Bayesian framework for calibration and credibility assessment of a 
computer simulator. The framework is particularly developed for simulators with 
uncertain subsystem models that are represented using functions. Figure [TJ shows 
schematic of the proposed framework for a complex system consisting of physically 
or mechanically interconnected subsystems. The system is investigated using ex- 
perimental observations and computer simulators, that use available information 
about the system for initial setup, while, repeated runs of the experiment and the 
simulator are used to understand more about the system. Experimental observa- 
tions of the system response are used in the Bayesian framework for calibration of 
the computer simulator. 

Most appealing facet of the Bayesian framework is its ability to provide the com- 
plete posterior statistics. However barring very simple cases, statistical sampling 
techniques are required for solution of the Bayesian calibration and uncertainty 
propagation problems. Markov Chain Monte Carlo (MCMC) method [13, G3] is 
one of the most widely used sampling technique for the Bayesian calibration. How- 
ever, exploration of the posterior distribution using the MCMC requires collection 
of a large number of samples for satisfactory approximation (often in the range 
of 10 3 — 10 6 ), rendering the Bayesian framework computationally prohibitive for 
a large scale system simulator. Thus, for generalized application of a Bayesian 
framework to the complex large scale system simulators, development of a compu- 
tationally efficient Bayesian calibration technique is essential. 

Marzouk et al. [n| have proposed a spectral projection based method for com- 
putationally efficient Bayesian calibration. The method uses a spectral expansion 
of a prior in generalized Polynomial Chaos basis. The Polynomial Chaos based 
spectral projection method is extensively investigated in the literature as a com- 
putationally efficient alternative to statistical methods for uncertainty propagation 
with comparable accuracy Polynomial chaos method is based on a concept 
of Homogeneous Chaos introduced by Wiener 21, where a random variable 
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is spectrally expanded in terms of Hermite polynomials. Cameron and Martin j23l ] 
have shown that any non-linear functional can be expanded in terms of a series 
of Hermite polynomials in L 2 sense. Although earlier attempts at using the poly- 
nomial chaos method (especially for turbulent fluid flow modeling) were not very 
successf ul |24M26ll , the method is found to be useful for solution of stochastic finite 
element |27l429l| and stochastic fluid flow problems [13, EiH . Xiu and Karniadakis 
32 . [33| ] have generalized the polynomial chaos method for spectral projection in 
terms of the Askey scheme of polynomials (HI]. The generalized Polynomial Chaos 
(gPC) method have been applied by various researchers for uncertainty propagation 
through simulators of systems of engineering importance [35- 38] . 



The method proposed by Marzouk et al. [19[ uses the gPC for propagation of 
the prior uncertainty to the simulator predictions. The resultant gPC expansion 
of the simulator predictions is used to define the likelihood. On availability of the 
experimental observations, probability distribution of the gPC basis is updated 
using the Bayesian calibration. On substitution of respective polynomial chaos 
coefficients, posterior distribution of parameters is obtained. Marzouk et al. [3^] 
have further extended the method for inference of spatially/temporally varying 
uncertain parameters. 

In this paper, the method proposed by Marzouk et al. [3^] is extended for the 
prior with uncertain hyper-parameters. Though a family of probability distribution 
to represent the prior uncertainty can be specified, associated hyper-parameters 
are rarely known deterministically. Realistic quantification of the prior uncertainty 
requires specification of probability distribution for uncertain hyper-parameters. 
Use of the uncertain hyper-parameters is more ubiquitous in case of calibration 
of simulators with model structural uncertainty. Hierarchical Bayesian inference is 
pro posed in the literature for calibration in presence of uncertain hyper-parameters 
[13, EH- However, methodology proposed by Marzouk et al. does not ex- 

plicitly consider the effect of uncertain hyper-parameters in the formulation. To 
make the spectral stochastic projection based Bayesian inference more precise, it 
is necessary to include uncertain hyper-parameters in the formulation. 



This paper proposes an extension of the method of Marzouk et al. [19|,|39[ to take 
into consideration the uncertainty in hyper-parameters of the prior distribution. A 
methodology is proposed to obtain Karhunnen-Loeve expansion (KL expansion) of 
a stochastic process in terms of functions of the hyper-parameters. The prior un- 
certainty in hyper-parameters is expanded in the gPC basis. Galerkin projection is 
used to evaluate gPC coefficients of the resultant KL expansion terms of a stochas- 
tic process. The prior uncertainty in subsystem model, represented in terms of 
gPC basis, is propagated to the simulator predictions using the intrusive Galerkin 
projection approach [29]. The Bayesian calibration is reformulated as a MCMC 
sampling from the posterior distribution of the gPC basis. The resultant gPC ex- 
pansion with posterior distribution of the basis defines the posterior distribution 
of the uncertain parameters and the model structure. The posterior distribution of 
the model structure defines credibility of the simulator model. The posterior pa- 
rameters are identified that quantifies acceptability of the simulator. The proposed 
method is demonstrated using a simulator of a quasi-one-dimensional flow through 
a nozzle. The particular choice of the application is motivated by the fact that the 
quasi-one-dimensional nozzle flow is well understood and can be simulated with 
limited computational resources. 

This research extends existing state of the art by: (a) extending the gPC expan- 
sion based method of Marzouk et al. (author?) [19i . |39J for priors with uncertain 
hyper-parameters; and (b) providing guidelines for acceptability of the simulator 
model using hyper-parameters of the posterior distribution. Note that a prelimi- 
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nary version of this work was reported in [4(| , while this article is significantly ex- 
panded by including (a) model structural uncertainty; (b) substantially elaborate 
theoretical analysis; and (c) additional numerical results to establish computational 
efficiency and efficacy of the method. 

The rest of the paper is organized as follows. Section 2 provides statistical for- 
mulation of the problem. In section 3, proposed method is discussed in detail. In 
section 4, numerical results for the calibration of a quasi-one-dimensional nozzle 
flow simulator are presented and finally in section 5, the paper is summarized and 
concluded. 



2. Statistical Formulation 

Let the system be investigated by observing M system responses, while, 
Tj(x, u(x s )) be an available simulator of the j th system response, where x G X 
is a set of deterministic control inputs and u(x s ) are uncertain subsystem model. 
Note that u(x s ) are functions with argument x s consisting of mix of some of the 
elements of x, some intermediate calculations from the system model and indepen- 
dent subsystem specific input parameters. Output of the subsystem model is used 
in the system model for further calculations. For brevity, the discussion presented 
in this paper assumes a single uncertain subsystem model; however, the proposed 
method can be extended to a more generic case without any change. Let Ui(x s ) 
denote the 'true' subsystem model, while, C/'( x > ut(x s )) be the 'true' but unknown 
j th system response prediction. Note that the simulator Tj(-, •) approximates the 
system response within limits of available knowledge. Thus conditional on the true 
subsystem model ut(x s ), Tj(x, iti(x s )) deviates from £.,-(x, u t (x s )) by 

£j(x,ut(x s )) = Tj(x,ut(-K s )) + <5j(x), (1) 

where <5j(x) is a discrepancy function. 

Let y e (x) represent an experimental observation of the system at a control input 
setting x, while 6j(x) denote the corresponding measurement uncertainty. Relation- 
ship between experimental observation and the true system response is given by 

y ej ( x ) = Cj(x,n t (x s )) + ej(x). (2) 

Let experimental observations are obtained at iV input conditions and denote a 
set of experimental observations by Y e = {Y e .; j = 1,...,M} where Y e . = 
{y ej (x.i); i = 1,...,N}. Similarly define 5j = {5j(xj); i = 1,...,N} and S = 
{Sj; j = 1, ...,M}. Also define u t = {u t (x Sti ); i = 1,...,N U }, where typically N u 
is significantly greater than N. Note that represent a realization of a random 
function Uf(pc s ) at N u input settings. Information available from the experimental 
observations is used in the Bayes theorem as 

V(u t , S | Y e ) oc V(Y e | u t , 6) x V(u t , 6), (3) 

where V(u t , S) is a prior, V(Y e \ u t , S) is likelihood and V(u t , S \ Y e ) is a posterior 
probability distribution. 

Prior probability in and S is specified using independent Gaussian processes. 
However, complete definition of the prior requires specification of uncertain hyper- 
parameters of the probability distribution. Let 9 U € @ u and 0$ € 0,5 are the 
hyper-parameters of subsystem model and discrepancy function respectively, where 
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® u and 05 are set of possible values. Thus in the presence of the uncertain hyper- 
parameters, posterior probability distribution ([3]) takes the form 

V(u t (9 u ), 6(0 S ), On, Os | Y e ) oc T(Y e \ u t (6 u ), 8{0 6 ), 9 U , S ) 

( 4 ) 

x v(u t (6 u ),d(0s) | e u ,e 5 ) x v{e u ) x v(e 5 ). 

For better readability, dependence of the uncertain subsystem model and the dis- 
crepancy function on uncertain hyper-parameters is explicitly shown henceforth. 

In the present paper, Bayesian inference is developed assuming the Gaussian 
process prior for discrepancy functions, as it is extensively used in the literature 
for specification of prior on random functions 1(3, 41-43)]. Prior uncertainty in sub- 



system model and discrepancy function is assumed to be independent. Uncertainty 
in the experimental observations is specified by a zero mean normally distributed 
random variable 6j with standard deviation a e . . The experimental uncertainty at 
different input settings x is assumed to be uncorrelated. On marginalization of 
Sj(x), posterior distribution is given by 

m i r i i 
v(u t (e u ),s(e u ),e u ,e 5 | Y e ) oc J]_== ea; p|-- (y e . -^yvr 1 (Y e . 



4 y/TZi 

x v(u t (e u ) | e u ) x v(p u ) x v{o s ), 



(5) 



where /J,j = {Tj(xi,u t (x s ;6 u )) + E(5j(xi;9 s )); i = 1, ...,iV} while Ej = Y, s . +ap N , 
Ijv being N x N identity matrix. 

Metropolis-Hastings algorithm [3, SH] is used to sample from the posterior dis- 
tribution ([5]). For each sample, evaluation of Tj(x,ut(x s ;9 u )), E" 1 and | T,j \ im- 
pose significant computational expenses on the Bayesian framework. In the present 
paper, a gPC expansion based method is proposed for computationally efficient 
evaluations of Tj(x, tit(x s ; U )), E" 1 and | Ej |. 



3. Generalized Polynomial Chaos Expansion of a Gaussian Process with 
Uncertain Hyper-parameters 

For brevity, the proposed spectral formulation is described first for a susbsystem 
model, u(x; U ), which is subsequently extended for a discrepancy functional The 
formulation is derived for a zero mean Gaussian process, though, it can easily be 
extended for non-zero mean processes. The derivation uses KL expansion of the 
Gaussian process with uncertain hyper-parameters, which is projected on a gPC 
basis using the intrusive Galerkin projection. 



3.1. KL Expansion 

Let u(x;9 u ) be a zero-mean Gaussian process with a covariance function 
C u (xi, X2; 6 U ). The covariance function can be approximated as |29|] 



N 

C u (x 1 ,x 2 ;O u ) = A n (fl„)e ra (xi, M )e n (x 2 , 6 U ), 

n=l 



(6) 



For notational convenience, x s is replaced by x in this section. 
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where N is the number of expansion terms retained in the spectral approximation. 
A n (#u) and e n (x, 6 U ) are eigenvalues and eigenfunctions of the covariance kernel, 
which are given by solution of the Fredholm's integral equation of the second kind 



/ C u (x 1 ,x- 2 ;6 u )e n (xi,6 u )dx 1 = X n (9 u )e n (x 2 ,9 u )- (7) 
Jx 

Explicit dependence of the eigenvalues and eigenfunctions on the hyper-parameters 
6 U should be noted. The resultant KL expansion using ([6]) is given by 29J] 



N 

u(x; U ) = ^2 \/%Jduje n (x, 6 u )x n , (8) 

n=l 



where Xn are independent zero-mean standard normal random variables. For a 
given 6 U , the eigenvalue problem Q can be numerically solved using a Galerkin 



projection based approach [47|]. In this paper, the approach is extended for uncer- 



tain hyper-parameters as follows. 
Eigenfunctions e n (x,6 u ) can be spectrally approximated as 



v 



e n (x,0 u ) = £< l (0 u )Vi(x), (9) 



i=l 



where ipi(x) are Legendre polynomials and d?(0 u ) are respective expansion coef- 
ficients. Use in ([7]), multiply both sides by tpj(x.2) and integrate w.r.t. dx.2 to 
obtain 



y^d?(6 u ) / / C u (xi,x 2 ;0 u )'i/'i(xi)'0y(x2)dxidx2 
~1 Jx Jx 

( 10 ) 

= A n (0 w )5^(0„) / V'i(x 2 )V'i(x2)rfx2. 



i=l Jx 



Using 

Aj(0 u ) = / C(xi,x 2 ; 0„)V'i(xi)V'j(x2)dxidx2, D ij {9 u ) = d{(9 u ), 
Jx Jx 

Bij{9 u ) = / Vi(x2)^j(x2)dx 2 , Ki(Qu) = A;(0 U ), 

(|10p can be written in a matrix form as 

A(0 u )D(0 u ) = A(0 U )B(6 U )D(0 U ), (11) 



which is a generalized eigenvalue problem (GEP) that can be solved using the QZ 
algorithm [291 ] . 

The matrices in (jlip are functions of 8 U . Thus, to solve (jlip . consider spectral 
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expansion of eigenvalues and eigenfunctions as 



A*(0«) = E l iU0u), dl(6 u ) = Ci,kH0u), (12) 



i=l i=l 



where (f>i(0 u ) are appropriately scaled Legendre polynomials that form complete 
orthonormal basis on L 2 (& u ). The coefficients are given by 



u J @ jUOu)do u > Cl ' k ~ ) & j 2 (o u )de u ■ 



Using the Gauss-Legendre quadrature, the integrals can be approximated as 

\ % {e u )U0u)de u = y^\ i {ei) ( j> i {ei)w q , 

(14) 

N„ 



d n k (e u )M0u)de u = Y J dl{ei)^{oi)^ 

9=1 



where N q are the number of quadrature points used, Q q u are the quadrature nodes 
while w q are the respective quadrature weights. The expansion coefficients, If and 
c" fe , are calculated by solving (fTT|) at quadrature nodes Q q u and substituting the 
solution in (|14p to evaluate integrals in (|f 3[) . 



3.2. ffPC Expansion of a Gaussian Process 
Hyper-parameters # u can be expanded in gPC basis as [3 

p 

e u = J2%H p (t), (is) 

where H p (£) are the Hermite polynomials, 6™ are respective expansion coefficients, 
while £ is a vector of independent standard normal random variables. Using the 
intrusive approach, gPC expansion of 4>i{9 u ) is given by 

p 

<p i {0u) = Y J ^v H v^)- (16) 

P =i 



Using (fT2|) and (fT6j) . gPC expansion of y \ n (Q u ) is given by 



L p 

Vm^) = E E (it) 
i=i p=i 



where s£ = (\/Z)i=i ^M^u), 4>k{0 u )) / (<f>l(9u)) , and (•, •) denote an inner prod- 
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uctd. Similarly, the gPC expansion of eigenfunctions can be obtained using ([12 
and ([TBI) as 



N L P 

e n (x, U ) = £££ c&^flpCO^Cx). (18) 
fe=i i=i P =i 

Thus, a zero-mean Gaussian process u(x; U ) can be expanded in gPC basis as 



u(x;fl tt ) = 5^«p(x)flp(0, (19) 

p=i 

where the expansion coefficients Ufc(x) are given by using (jHJ), (ITTj) and (fT8l) as 

n=lm=l i=l j=lp=l g=l \ k\*J/ 

(20) 

Note that since it(x; U ) is a Gaussian process, Xn in © are standard normal 
variables, thus, \ n = i/ n+ i(£) is used in ([201) . 



4. Stochastic Spectral Projection based Bayesian Calibration 

The prior uncertainty in the subsystem model, ut(x s ;6 u ), and the discrepancy 
function, <5j(x; 0$), is given by independent Gaussian processes. Using (fl9l) . spectral 
expansion of the prior is given by 

P P 
u t (x s ;O u ) =J2ti P (xs)H p (£); SjfrOs) = J^x)^). (21) 

p=i p=i 

Intrusive Galerkin projection approach is used to propagate the prior uncertainty 
in ut(x s ;6 u ) to the simulator predictions, thus, 

p 

TfautfaO*)) =£z?(x)fl p (0. (22) 
p=i 

Using the gPC expansion of 2}(x, u t (x s ; U )) and 5j(x; 6,5) in ([T]) to obtain the gPC 
expansion of (j(x,u t (x s ; U )) as 

p 

C,-(x,Ut(x s ;0„)) = 2gx)ff p ({), (23) 
p=i 



2 For example, inner product of functions /(■) and g(-) is given by 

</(-),9(-)>= / /(09(0**(0. 

where ^t(-) is a measure on L 2 (0) 
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where Cp( x ) = ^p( x ) + 5p(x). Use the gPC expansion ([23]) in © to obtain 

p 

y e3 (x) = ^g(x)if p (0 + e r (24) 
p=i 



Note that £ are the only uncertain variables in (|24p . thus, the Bayesian inference 
problem is reformulated as sampling from the posterior distribution of 

Let t = {3jj(x); p = 1,...,P; j = 1,...,M} and 5 = {$(x); p = 1,...,P; j = 
1, M} define the sets of respective gPC coefficients. Using the gPC expansion of 
the uncertain variables ( (|21J) - (124[) ) in (J^J), the Bayesian calibration problem can 
be reformulated in terms of £ that capture all the randomness in the system model 
and hyper-parameters: 

V{£ | Y e , f , 8) cx V(Y e | e, T, 5) x (25) 

Since £ are i.i.d. standard normal random variables, with independent Gaussian 
uncertainties in experimental data, the posterior distribution of £ can be obtained 
by: 

M N d 

P(£ I Y e , t, 6) cx JJ -= exp {-I (Y £j - H f EJ 1 (Y Cj - H ) } x [j e ~«/ 2 , (26) 



j=l V ' n=l 

where /x, = {£j =1 ^'( x i)# P (£) + ^(x 4 ); i = 1, • • • , N}. 

Markov Chain Monte Carlo (MCMC) method is used to sample from (|26|) . Note 
that MCMC applied to (I26p does not require solution of the simulation model, 
Tj(x,u(x s ,6 u )), thus, the posterior distribution can be explored efficiently. How- 
ever, solution of (|26p requires numerical evaluation of |Ej| and EJ 1 , which may 
impose non-trivial computational cost on the MCMC sampling. In the present 
paper, numerical evaluation of |Ey| and EJ 1 is accelerated as follows. 

The numerical evaluation of IS J and E7 1 is accelerated using the gPC expansion 
of the individual elements of the covariance matrix Ej(xi,X2) as 

p 

E i (x 1 ,x 2 ) = 5^qj(xi > x 2 )fl p (0. (27) 
p=i 

The determinant |S 7 -| and the individual elements of the inverse E" 1 can also be 
expanded in gPC basis as 

p p 
|E,-| =Y,D J p H p (£), Sj^xi.xa) = ^ ^( Xl ,x 2 )^(^), (28) 
p=i p=i 

where the gPC expansion coefficients are given by 

1-1 
"(Hi) 



Gaussian quadrature is used to evaluate polynomial chaos coefficients. (E ■ , 

and (|Ej|, Hk(£)) are calculated by evaluating E" 1 and |Ej| at quadrature nodes, 
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while Ii and Di are evaluated using (|29p . The resultant gPC expansion (|28p is 
used in ([26]) for MCMC sampling. 



5. Numerical Example: Calibration of Quasi-One-Dimensional Nozzle Flow 
Simulator 

5.1. Problem Setup 

A quasi-one-dimensional supersonic flow through a nozzle is considered to verify the 
proposed Bayesian calibration method. The quasi-one-dimensional flow is uniform 
across the cross-section with properties varying in x-direction, while, the effect 
of change in the cross sectional area is considered. The flow is defined using the 
compressible Euler equations in conservative form 

dq di 



m + d- x =^ (30) 



where, 




pvA \ 

q = | pvA | , f = | P v 2 A + PA , g=(p| 

pvEA + PvA) \ 

Here, p is density, v is velocity, A is cross sectional area, P is static pressure and 
E is the total energy per unit mass. In the governing equations, static pressure is 
substituted by the total energy per unit mass using 

£ =(^+r 2 ' < 31 > 

while the ideal gas equation is used for the closure. 

Variation of nozzle area A is assumed to be uncertain, which is inferred using 
the Bayesian framework. Stationary Gaussian process with known mean profile is 
used as a prior for uncertain nozzle area. The procedure of obtaining the stochastic 
spectral formulation is summarized as follows. Define q\ = pA, qi = pvA, q% = pEA 
and the corresponding gPC expansions 



p=i p=i p=i 

Multiplying both the sides by Hk($) and taking the inner product, governing equa- 
tions of a stochastic quasi-one-dimensional nozzle flow is given by 

where 



(Q), 
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7(7-1) 
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EP v^-P * <- 

i=l L,j=i l^k=i Q2,iQ2,jr k 



7-1 

7 



(HiH jH k H p ) . 

w> + 

(Hj Hj H k Hp ) 
{HI) 



Q3, P - Ei=i Ej=i Efc=i Q2,m,jr k 

V JV P V 7V P * * ~ {HiHjH k Hp) 
7 Z^i=i 2^=1 Z^fc=i Q2,iQ3,jr k ^ 



11 



^P ^P ^P v^P v^P 

H=l 2-^=1 l^k=l 1^1=1 l^m=l Q2,iq2,jq2,krir m 



(Hj Hj Hk Hi H m Hp ) 
(HI) 



\ 



7(7-1) 

27 



ly^P T P 



7 

=iE 



=1 Efc= 



1 Ef=i hthjhdA i dx &&$$EL 



V / 

and r = Note that the governing equations is a set of 3 x P partial differential 
equations. These governing equations are numerically solved using the central dif- 
ference scheme in the spatial dimension and the fourth order Runge-Kutta method 
in the temporal dimension. A uniform grid with Ax = 0.01 is used in the spatial 
dimension, while, the time step At = 0.0001 is used for time integration. The in- 
ner products involved in the governing equations are evaluated a-priory using the 
Gauss-Hermite quadrature. 

For demonstration purpose, the proposed method is applied using 'hypothetical 
test bed nozzle' data. The 'hypothetical nozzle' can be created using the nozzle 
simulator with sensors that can measure the nozzle response with typical accuracies. 
A hypothetical nozzle is specified through a set of parameters and 'true' area profile, 
all of which are treated as precisely known. In all the test cases presented in this 
section, density p, velocity v, pressure P and static temperature T are used as 
responses of interest. Steady state predictions using known nozzle area profile are 
used as experimental observations. The experimental uncertainty is specified by a 
zero-mean normally distributed random variable with standard deviation given by 
1% of the mean value. For all the test cases, inflow conditions are Mach number 
Mi n = 1.5, static pressure Pj n = 1.0 and density p = 1.0. Figure [2] shows a 
typical nozzle area profile and spatial variation of normalized responses at steady 
state. Note that supersonic flow through a divergent nozzle results in increased 
velocity, while the static pressure, temperature and density are decreased. As can 
be observed from Figure El this behavior is captured well by the nozzle simulator. 



5.2. Prior Uncertainty Propagation 

Prior uncertainty in the nozzle area is specified using a Gaussian process with 
known mean and the covariance function is given by 

C(xi,x 2 ) = <r 2 exp (—X(xi -x 2 ) 2 ) , (33) 

where a 1 is the variance and A is the correlation length, a 2 and A are treated as 
uncertain. Prior uncertainty in a 1 is specified using an inverse Gamma distribu- 
tion, IG(9. 0,0.5), while Gamma distribution, G(5. 0,0.2), is used as a prior for A. 
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Figure 2. Figure shows (a) Nozzle area variation and (b) steady state response predictions for deterministic 
case. All the steady state responses are normalized using inflow values. 
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Figure 3. Figure shows a) Comparison of prior mean nozzle area with true nozzle area and b) prior 
probability distribution of nozzle area. 



Figure [3] a) shows prior mean nozzle area while Figure E^b) shows prior probability 
distribution. Comparison of prior mean with 'true' nozzle area is also shown in 
Figure El[a) . 

The Bayesian framework is implemented using the intrusive gPC expansion for 
propagation of the prior uncertainty to the system response. The prior uncertainty 
is projected on a Hermite polynomial chaos basis. To investigate the trade-off be- 
tween computational efficiency and accuracy, results of polynomial chaos method 
for uncertainty propagation are compared with the Monte Carlo simulation. Total 
10000 samples are used for the Monte Carlo method. All the computations are per- 
formed on a desktop computer with intel i5-460 processor. CPU time is estimated 
using a FORTRAN intrinsic cpu_time routine. Figure H] shows the comparison of 
computational time requirement as a function of number of eigenfunctions used, 
N, for different polynomial chaos order, p. As can be observed from the figure, the 
computational cost increases polynomially with N. In particular, for a system with 
stochasticity of order s, the computational cost increases as Ns pq , where q is the 
order of non-linearity of the system and p is order of the gPC basis. 
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Figure 4. Figure shows cpu time required as a function of number of eigenfunctions used. Effect of the 
polynomial chaos order is also shown. 
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Figure 5. Figure shows (a) Li-error in mean and (b) Li-error in variance. Effect of the polynomial chaos 
order is also shown. 

Figure [5] shows the L\ -error in the mean and the variance. The L\ -error is cal- 
culated with respect to the Monte Carlo method. With the increase in number 
of eigenfunctions used and the polynomial order, error in the mean and the vari- 
ance reduces. For N = 4 and the 2 nd order polynomial chaos, method provides 
prediction of system response with error of the order of 10 -3 in both mean and 
variance at 10-times lower computational cost. Note that the computational cost 
of the proposed Bayesian calibration method is dominated by the computational 
time requirement for solution of forward propagation problem, while the computa- 
tion cost of the MCMC sampling is negligible. Thus, the conclusions drawn from 
the computational cost comparison for forward propagation can be extended to 
solution of the inverse problem without any change. 



5.3. Spectral Projection- Based Bayesian Calibration 

5.3.1. Baseline Case 

The computational cost of the proposed method is compared with the direct 
MCMC sampling for the Bayesian calibration. Prior for the uncertain nozzle area 
is specified using the Gaussian process with the uncertain variance and the co- 
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variance length. 7(7(9.0, 0.5) prior is used for the variance while G(5.0, 0.2) prior 
is used for the correlation length. First four eigenfunctions are used in the spec- 
tral expansion, while 2 nd order Hermite polynomials are used as the gPC basis. 
Model structural uncertainty is quantified by IG(9.0, 0.5) prior for the variance 
and G(6.0, 2.0) prior for the correlation length of the covariance function of the 
discrepancy function. The gPC expansion coefficients of the simulator output, ob- 
tained using the intrusive Galerkin projection, are used in the likelihood function 
to define the posterior distribution, which is explored using the MCMC sampling. 
Gaussian distribution centered on the present state is used as a proposal distri- 
bution for the Markov Chain. Total 100000 samples are collected after rejecting 
the initial 10000 samples. Results of the proposed method are compared with the 
direct implementation of MCMC for the Bayesian calibration, where the simulator 
output is used to define the likelihood. Total 10000 samples are collected using 
the direct MCMC after the initial burnout period of 1000 samples. The total CPU 
time required for the implementation of the Bayesian framework using the direct 
MCMC sampling is 12732.65 seconds. The CPU time requirement for the pro- 
posed gPC expansion based Bayesian framework is split into the time required for 
the intrusive Galerkin projection of the prior uncertainty to the system response 
and the MCMC sampling from the posterior distribution. For the present test 
case, intrusive Galerkin projection is implemented in 1194.15 seconds, while the 
total CPU time for the MCMC sampling is 8.48 seconds, taking 1202.63 seconds 
for complete implementation of the proposed stochastic spectral projection based 
Bayesian framework. 

Figure [6] shows comparison of the posterior mean nozzle area obtained using the 
direct MCMC and the proposed method. The prior mean nozzle area is also shown 
in the figure. The posterior mean nozzle area obtained using both the methods 
matches closely with the 'true' nozzle area. The comparison of the posterior dis- 
tribution for the hyper-parameters of the uncertain nozzle area is shown in Figure 
[3 The probability distribution of the hyper-parameters is not updated noticeably 
after the Bayesian calibration, as the experimental observations of the system re- 
sponse are not expected to contain significant information about the covariance 
structure of the uncertain nozzle area. From the figures, it may be concluded that 
the proposed spectral projection based Bayesian framework provide inference of 
the uncertain parameters with accuracy comparable to the implementation of the 
direct MCMC sampling at significantly lower computational cost. 

5.3.2. Choice of Discrepancy Function 

To investigate the effect of choice of the prior for discrepancy function, Bayesian 
framework is implemented using different priors for variance. Prior uncertain in 
the nozzle area is specified as discussed earlier. Prior uncertainty in discrepancy 
function is specified using zero-mean Gaussian process with squared exponential 
covariance function (|33|) . where the variance and the correlation length are un- 
certain hyper-parameters. For all the test cases presented in this section, prior 
in correlation length is represented using G(6.0, 2.0), specifying correlated discrep- 
ancy function. For variance, IG(6.0, 2.0) and 7(7(1.5, 2.0) priors are investigated. To 
simulate model structure uncertainty, hypothetical test bed data is obtained using 
viscous model, whereas, simulator is defined using inviscid model. Total 5 testbed 
data points are used with 1% standard deviation. Figure [8] shows comparison of pos- 
terior mean nozzle area with true nozzle area and the prior mean. Posterior mean 
for 1(7(6.0, 2.0) prior matches closely with the true nozzle area, whereas, posterior 
mean nozzle area for 7(7(1.5, 2.0) prior deviates from the true nozzle area. Figure 
[8] demonstrates the significant impact of the prior model structural uncertainty on 
the inference of the uncertain parameters. Confidence on the available model is 
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Figure 6. Comparison of the posterior mean nozzle area with the true nozzle area and prior mean. The 
comparison is shown for the posterior mean nozzle area obtained using implementation of the direct MCMC 
sampling and the proposed method. 
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Figure 7. Figure shows comparison of posterior distribution for the hyper-parameters of the uncertain 
nozzle area obtained using the direct MCMC and the proposed spectral projection based method. Figure 
(a) shows the comparison for the variance and the Figure (b) shows the comparison for the correlation 
length. 

specified through the prior on the variance. In the case of high confidence on the 
simulator model, specified through the high value of a for IG(a, (3) prior, significant 
amount of the information provided by the data is used to update the parameter 
uncertainty. If the low confidence prior is specified for the model through the lower 
value of a for IG(a, /3) prior, the calibration process uses significant information 
provided by the data to improve confidence on the simulator model, whereas, less 
information is used for updating the uncertain parameters. Note that in the case 
of very high confident prior on the model structure (approaching the scenario of 
no model structural uncertainty), calibration process attributes any remnant error 
in the model to the parameters. Thus, the present authors propose to avoid priors 
that specify either low confidence or very high confidence on the model structure. 
In the remaining test cases presented in this paper, /G(6.0, 2.0) prior is used for 
the variance of the discrepancy function. 



5. 3. 3. Effect of Simulator Model Error 



To investigate efficacy of the calibration process in presence of error in simulator 
model, artificial discrepancy is introduced in the simulator model by multiplying 
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Figure 8. Figure shows effect of the choice of prior for discrepancy function on posterior mean of the 
nozzle area. Note that 7G(1.5, 2.0) represent lower confidence on the simulator model as compared to the 
JG(6.0,2.0) prior. 

the gPC expansion coefficients of the static pressure by f .5. Modified chaos coeffi- 
cients are used in the Bayesian framework. The framework is implemented using fol- 
lowing two test cases: (1) simulator model without taking into account discrepancy 
function (which signifies full confidence on the simulator model), and (2) G(6.0, 2.0) 
and 1(7(6.0, 2.0) priors for correlation length and variance of the discrepancy func- 
tion (\$ and cr| respectively). Figure [9] shows comparison of posterior nozzle area 
with the true nozzle area. For calibration without considering discrepancy func- 
tion, Bayesian framework assumes simulator to be the 'true' representation of the 
physics. The Bayesian calibration method attributes all the difference between the 
test bed data and the simulator prediction to the uncertain parameters. Thus, re- 
sultant posterior mean deviates significantly from the true nozzle area. However, 
when discrepancy function is considered in the Bayesian calibration, no significant 
update is observed in the nozzle area. 

Figure [10] shows comparison of the prior and the posterior probability distri- 
bution of a\ and A,5. The posterior probability distribution of 1 / cr| for the static 
pressure is shifted towards left, indicating the lower posterior expected value of 
1/cj| as compared to the prior. However, no significant change is observed for the 
probability distribution of A5. Note that the prior G(6.0, 2.0) indicated correlated 
discrepancy in the simulator model. Thus, the posterior probability distribution 
of the discrepancy function for the static pressure indicates correlated discrepancy 
in the simulator model with the high expected value of <r|. As the error in the 
experimental observations enters the Bayesian framework as the uncorrelated un- 
certainty, the highly correlated discrepancy is expected to result from the error 
in the simulator model. Thus, the posterior distribution of the discrepancy func- 
tion signals the need for verification and validation of the simulator, particularly 
subroutines affecting the static pressure. 

5.3.4- Effect of Erroneous Observations 

To investigate effect of the erroneous experimental observations on the Bayesian 
framework, proposed method is implemented with artificially introduced discrep- 
ancy in test bed data. The methodology is demonstrated by multiplying 1.5 to 
static pressure data. Figure [TT] shows resultant posterior mean nozzle area. As ob- 
served in case of artificially introduced discrepancy in the simulator, significant 
deviation is observed between posterior mean and true nozzle area when discrep- 
ancy function is not considered in the formulation. However, when discrepancy 
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Figure 9. Figure shows the effect of error in the simulator model on the posterior nozzle area. Results 
are obtained for artificially introduced discrepancy in the simulator model. The comparison is shown for 
the Bayesian calibration without using the discrepancy function (no discr) and with using the discrepancy 
function (w/ discr). 
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Figure 10. Figure shows comparison of prior and posterior distribution for a) and b) \g in presence of 
error in the simulator model. Results are obtained for artificially introduced discrepancy in the simulator 
model. 



function is considered, no significant change is observed in the nozzle area. 

Figure [12] shows posterior probability distribution for a\ and \$. Posterior dis- 
tribution of 1/cr? for static pressure has moved towards left, indicating the high 
posterior discrepancy. Posterior probability distribution of As for static pressure 
is shifted towards right, indicating highly uncorrelated discrepancy. Since the un- 
certainty due to unknown or poorly known physics is expected to result in the 
correlated discrepancy, uncorrelated discrepancy may be attributed to the experi- 
mental observations. Thus the posterior distribution indicates need for the review 
of the experimental observations. 

5. 3. 5. Summary 

The numerical test cases presented in this paper have demonstrated ability of 
the proposed Bayesian framework to infer spatially /temporally varying uncertain 
parameters in the presence of model structural uncertainty. In addition to the in- 
ference of the uncertain parameters, the proposed Bayesian framework provides 
useful insight into the credibility of the simulator model. The posterior probability 
distributions of u? and A,5, through hyper-parameters a and /3, are identified as in- 
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Figure 11. Figure shows the effect of erroneous experimental observations on the posterior nozzle area. 
Results arc obtained for artificially introduced errors in the experimental observations. 
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Figure 12. Figure shows comparison of prior and posterior distribution for a) <r| and b) in presence of 
erroneous experimental observations. 



dicators of the veracity and validity of the simulator model. Based on the posterior 
values of the hyper-parameters a and /3 for the posterior probability distributions 
(a CT5 and f3 as ) and A,5 (a\ s and authors provide following guidelines: 



(1) If the posterior a Us is greater than the prior a Us , then the simulator model 
should be accepted with improved confidence. 

(2) a as < 1 : The posterior distribution does not have mode resulting in maxi- 
mum probability for \ — > 0. For such posterior, calibrated simulator model 
should be rejected. 

• If mean and mode of posterior distribution of A$ indicate strong corre- 
lation for discrepancy function (typically a a < 1 or (3\ > a\), rigorous 
verification and validation process is advised with a focus on subsystem 
models that predict system responses for which a as < 1. 

• If mean and mode of posterior distribution of A$ indicate very weak 
correlation for discrepancy function (typically a\ > 1 or (3\ « a), 
review of experimental observations is advised. 

• For all the other cases, rigorous verification and validation of simulator 
model is advised with a note on review of experimental observations. 



(3) 



On 



> 1 and p as > 



Calibrated simulator can be used, however, high 
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uncertainty in prediction of the simulator response should be expected. As 
per discussion void point (2), verification and validation of simulator model 
and a review of experimental observations is advised. 
(4) a as >> 1 and (3 ag < a ag : Calibrated simulator can be used with high 
confidence. 



6. Concluding Remarks 

This paper has demonstrated computational efficiency of a gPC based Bayesian 
framework for calibration of a large scale system simulator. The proposed frame- 
work has extended the established method to the priors with uncertain hyper- 
parameters. Efficacy of the proposed Bayesian framework is demonstrated for cal- 
ibration of a quasi-one-dimensional divergent nozzle flow simulator. Ability of the 
method to infer spatially /temporally varying uncertain parameters is shown us- 
ing the update of the nozzle area. The proposed method has provided accurate 
inference of nozzle area at one-tenth of a computational cost as compared to the 
direct implementation of the Bayesian framework. Hyper-parameters of the pos- 
terior distribution of model structure uncertainty are identified that provide in- 
formation about veracity and validity of the computer simulator. Based on the 
hyper-parameters, guidelines have been provided for acceptability of the simulator 
model. Although demonstrated for a specific set of priors, the proposed method is 
generic in nature and can admit arbitrary priors. However, depending on the gPC 
basis used, higher order polynomials may be required for the satisfactory spectral 
approximation of the prior, incurring comparatively higher computational cost on 
the Bayesian framework. 
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